
function [Wff_h,Wfi_h,Wfn_h,Wif_h,Wnf_h,Wii_h,Win_h,Wni_h,Wnn_h,Wff_l,Wfi_l,Wfn_l,Wif_l,Wnf_l,Wii_l,Win_l,Wni_l,Wnn_l,iterVF]=solve_value_function_A1(Par1,Ff_1,wf_1,Fi_1,wi_1,Ff_2,wf_2,Fi_2,wi_2,n)

r=0.036;

% Clenshaw-Curtis quadrature
% n+1 nodes: x in [-1,1] by decreasing order
[x,wcc] = ccquad(n-1); x=flipud(x);

transition_rates=Par1(1:16)

% h: healthy, l: unhealthy
theta=Par1(21)
a_h=Par1(22)
gamma=0; % gamma = 0 at benchmark
b1_h = Par1(23); b2_h = Par1(24);
a_l=Par1(25)
b1_l = Par1(26); b2_l = Par1(27);
nu_hl= Par1(28); nu_lh = Par1(29);

df_1=transition_rates(1);
di_1=transition_rates(2);
lnf_1=transition_rates(3);
lni_1=transition_rates(4);
lfi_1=transition_rates(5);
lif_1=transition_rates(6);
df_2=transition_rates(7);
di_2=transition_rates(8);
lnf_2=transition_rates(9);
lni_2=transition_rates(10);
lfi_2=transition_rates(11);
lif_2=transition_rates(12);
p_1=transition_rates(13);
q_1=transition_rates(14);
p_2=transition_rates(15);
q_2=transition_rates(16);

% consider preallocation to speed up


%% Transforming the probabilities g and f in densities
% If wages are linearly spaced
wht_f_1=(wf_1(n)-wf_1(1))/(n-1);
wht_i_1=(wi_1(n)-wi_1(1))/(n-1);
wht_f_2=(wf_2(n)-wf_2(1))/(n-1);
wht_i_2=(wi_2(n)-wi_2(1))/(n-1);

% Initializing

% value function (min and max)
Wnn_h=(1-exp(-theta*(min(wi_1)+min(wi_2))))/theta/r;
Wnn_l=(1-exp(-theta*(min(wi_1)+min(wi_2))))/theta/r;
Wfn_min= (1-exp(-theta*(min(wf_1)+min(wi_2))))/theta/r; Wfn_max= (1-exp(-theta*(max(wf_1)+max(wi_2))))/theta/r;
Win_min=(1-exp(-theta*(min(wi_1)+min(wi_2))))/theta/r; Win_max=(1-exp(-theta*(max(wi_1)+max(wi_2))))/theta/r ;
Wnf_min=(1-exp(-theta*(min(wi_1)+min(wf_2))))/theta/r ; Wnf_max=(1-exp(-theta*(max(wi_1)+max(wf_2))))/theta/r ;
Wni_min=(1-exp(-theta*(min(wi_1)+min(wi_2))))/theta/r ; Wni_max=(1-exp(-theta*(max(wi_1)+max(wi_2))))/theta/r ;
Wff_min=(1-exp(-theta*(min(wf_1)+wf_2)))/theta/r; Wff_max= (1-exp(-theta*(max(wf_1)+wf_2)))/theta/r;
Wii_min=(1-exp(-theta*(min(wi_1)+wi_2)))/theta/r ; Wii_max= (1-exp(-theta*(max(wi_1)+wi_2)))/theta/r ;
Wfi_min=(1-exp(-theta*(min(wf_1)+wi_2)))/theta/r ; Wfi_max=(1-exp(-theta*(max(wf_1)+wi_2)))/theta/r;
Wif_min=(1-exp(-theta*(min(wi_1)+wf_2)))/theta/r ; Wif_max=(1-exp(-theta*(max(wi_1)+wf_2)))/theta/r;

Wfn_h=(Wfn_min+Wfn_max)/2+x*(Wfn_max-Wfn_min)/2; % Interpolation with CC quadrature (defined above)
Win_h=(Win_min+Win_max)/2+x*(Win_max-Win_min)/2;
Wnf_h=(Wnf_min+Wnf_max)/2+x*(Wnf_max-Wnf_min)/2;
Wni_h=(Wni_min+Wni_max)/2+x*(Wni_max-Wni_min)/2;
Wff_h=(Wff_min+Wff_max)*ones(1,n)/2+(Wff_max-Wff_min)*x'/2;
Wii_h=(Wii_min+Wii_max)*ones(1,n)/2+(Wii_max-Wii_min)*x'/2;
Wfi_h=(Wfi_min+Wfi_max)*ones(1,n)/2+(Wfi_max-Wfi_min)*x'/2;
Wif_h=(Wif_min+Wif_max)*ones(1,n)/2+(Wif_max-Wif_min)*x'/2;

Wfn_l=(Wfn_min+Wfn_max)/2+x*(Wfn_max-Wfn_min)/2; % Interpolation with CC quadrature (defined above)
Win_l=(Win_min+Win_max)/2+x*(Win_max-Win_min)/2;
Wnf_l=(Wnf_min+Wnf_max)/2+x*(Wnf_max-Wnf_min)/2;
Wni_l=(Wni_min+Wni_max)/2+x*(Wni_max-Wni_min)/2;
Wff_l=(Wff_min+Wff_max)*ones(1,n)/2+(Wff_max-Wff_min)*x'/2;
Wii_l=(Wii_min+Wii_max)*ones(1,n)/2+(Wii_max-Wii_min)*x'/2;
Wfi_l=(Wfi_min+Wfi_max)*ones(1,n)/2+(Wfi_max-Wfi_min)*x'/2;
Wif_l=(Wif_min+Wif_max)*ones(1,n)/2+(Wif_max-Wif_min)*x'/2;

critVF=1;
iterVF=1;
omegaW=0.1;
while mean(critVF)>0.005 && iterVF<=400
    
    % Reservation wages
    for i=1:n
        % to construct Wfn_h
        dWni_Wfn=abs(Wni_h-Wfn_h(i))./Wfn_h(i);
        Rw_ni_fn_h(i)=wi_2(min(find(dWni_Wfn<=min(dWni_Wfn))));
        Fi_2_Rw_ni_fn(i)=Fi_2(min(find(dWni_Wfn<=min(dWni_Wfn))));
        dWin_Wfn=abs(Win_h-Wfn_h(i))./Wfn_h(i);
        Rw_in_fn_h(i)=wi_1(min(find(dWin_Wfn<=min(dWin_Wfn))));
        Fi_1_Rw_in_fn(i)=Fi_1(min(find(dWin_Wfn<=min(dWin_Wfn))));
        dWff_Wfn=abs(Wff_h(i,:)-Wfn_h(i))./Wfn_h(i);
        Rw_ff_fn_h(i)=wf_2(find(dWff_Wfn<=min(dWff_Wfn), 1 ));
        Ff_2_Rw_ff_fn(i)=Ff_2(find(dWff_Wfn<=min(dWff_Wfn), 1 ));
        dWnf_Wfn=abs(Wnf_h-Wfn_h(i))./Wfn_h(i);
        Rw_nf_fn_h(i)=wf_2(find(dWnf_Wfn<=min(dWnf_Wfn), 1 ));
        Ff_2_Rw_nf_fn(i)=Ff_2(find(dWnf_Wfn<=min(dWnf_Wfn), 1 ));
        Ff_2_minRw_ff_fn_nf_fn(i)=Ff_2(min(find(dWnf_Wfn<=min(dWnf_Wfn), 1 ),min(find(dWff_Wfn<=min(dWff_Wfn)))));
        dWff_Wnf=abs(Wff_h(i,:)'-Wnf_h)./Wnf_h;
        Rw_ff_nf(i)=wf_2(find(dWff_Wnf<=min(dWff_Wnf), 1 ));
        Ff_2_Rw_ff_nf(i)=Ff_2(min(find(dWff_Wnf<=min(dWff_Wnf))));
        dWfi_Wfn=abs(Wfi_h(i,:)-Wfn_h(i))./Wfn_h(i);
        Rw_fi_fn(i)=wi_2(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
        Fi_2_Rw_fi_fn(i)=Fi_2(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
        
        derWfn_h(i)=1./(r+df_1*(1-p_2)+df_1*p_2*Fi_2_Rw_ni_fn(i)+lfi_1*Fi_1_Rw_in_fn(i) ...
            +lnf_2*Ff_2_minRw_ff_fn_nf_fn(i)+lni_2*Fi_2_Rw_fi_fn(i));
        
        dWni_Wfn=abs(Wni_l-Wfn_l(i))./Wfn_l(i);
        Rw_ni_fn_l(i)=wi_2(min(find(dWni_Wfn<=min(dWni_Wfn))));
        Fi_2_Rw_ni_fn(i)=Fi_2(min(find(dWni_Wfn<=min(dWni_Wfn))));
        dWin_Wfn=abs(Win_l-Wfn_l(i))./Wfn_l(i);
        Rw_in_fn_l(i)=wi_1(min(find(dWin_Wfn<=min(dWin_Wfn))));
        Fi_1_Rw_in_fn(i)=Fi_1(min(find(dWin_Wfn<=min(dWin_Wfn))));
        dWff_Wfn=abs(Wff_l(i,:)-Wfn_l(i))./Wfn_l(i);
        Rw_ff_fn_l(i)=wf_2(find(dWff_Wfn<=min(dWff_Wfn), 1 ));
        Ff_2_Rw_ff_fn(i)=Ff_2(find(dWff_Wfn<=min(dWff_Wfn), 1 ));
        dWnf_Wfn=abs(Wnf_l-Wfn_l(i))./Wfn_l(i);
        Rw_nf_fn_l(i)=wf_2(find(dWnf_Wfn<=min(dWnf_Wfn), 1 ));
        Ff_2_Rw_nf_fn(i)=Ff_2(find(dWnf_Wfn<=min(dWnf_Wfn), 1 ));
        Ff_2_minRw_ff_fn_nf_fn(i)=Ff_2(min(find(dWnf_Wfn<=min(dWnf_Wfn), 1 ),min(find(dWff_Wfn<=min(dWff_Wfn)))));
        dWff_Wnf=abs(Wff_l(i,:)'-Wnf_l)./Wnf_l;
        Rw_ff_nf(i)=wf_2(find(dWff_Wnf<=min(dWff_Wnf), 1 ));
        Ff_2_Rw_ff_nf(i)=Ff_2(min(find(dWff_Wnf<=min(dWff_Wnf))));
        dWfi_Wfn=abs(Wfi_l(i,:)-Wfn_l(i))./Wfn_l(i);
        Rw_fi_fn(i)=wi_2(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
        Fi_2_Rw_fi_fn(i)=Fi_2(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
        
        derWfn_l(i)=1./(r+df_1*(1-p_2)+df_1*p_2*Fi_2_Rw_ni_fn(i)+lfi_1*Fi_1_Rw_in_fn(i) ...
            +lnf_2*Ff_2_minRw_ff_fn_nf_fn(i)+lni_2*Fi_2_Rw_fi_fn(i));
        
        % to construct Win_h
        dWni_Win=abs(Wni_h(i,:)-Win_h(i))./Win_h(i);
        Rw_ni_in_h(i)=wi_2(min(find(dWni_Win<=min(dWni_Win))));
        Fi_2_Rw_ni_in(i)=Fi_2(find(dWni_Win<=min(dWni_Win), 1 ));
        dWfn_Win=abs(Wfn_h(i,:)-Win_h(i))./Win_h(i);
        Rw_fn_in_h(i)=wf_1(find(dWfn_Win<=min(dWfn_Win), 1 ));
        Ff_1_Rw_fn_in(i)=Ff_1(find(dWfn_Win<=min(dWfn_Win), 1 ));
        dWif_Win=abs(Wif_h(i,:)-Win_h(i))./Win_h(i);
        Rw_if_in_h(i)=wf_2(find(dWif_Win<=min(dWif_Win), 1 ));
        Ff_2_Rw_if_in(i)=Ff_2(find(dWif_Win<=min(dWif_Win), 1 ));
        dWnf_Win=abs(Wnf_h(i,:)-Win_h(i))./Win_h(i);
        Rw_nf_in_h(i)=wf_2(find(dWnf_Win<=min(dWnf_Win), 1 ));
        Ff_2_Rw_nf_in(i)=Ff_2(find(dWnf_Win<=min(dWnf_Win), 1 ));
        Ff_2_minRw_if_in_nf_in(i)=Ff_2(min(find(dWif_Win<=min(dWif_Win), 1 ),min(find(dWnf_Win<=min(dWnf_Win)))));
        dWif_Wnf=abs(Wif_h(i,:)'-Wnf_h)./Wnf_h;
        Rw_if_nf_h(i)=wf_2(min(find(dWif_Wnf<=min(dWif_Wnf))));
        Ff_2_Rw_if_nf(i)=Ff_2(min(find(dWif_Wnf<=min(dWif_Wnf))));
        dWii_Win=abs(Wii_h(i,:)-Win_h(i))./Win_h(i);
        Rw_ii_in_h(i)=wi_2(min(find(dWii_Win<=min(dWii_Win))));
        Fi_2_Rw_ii_in(i)=Fi_2(min(find(dWii_Win<=min(dWii_Win))));
        
        
        derWin_h(i)=1./(r+di_1*(1-q_2)+di_1*q_2*Fi_2_Rw_ni_in(i)+lif_1*Ff_1_Rw_fn_in(i) ...
            +lnf_2*Ff_2_minRw_if_in_nf_in(i)+lni_2*Fi_2_Rw_ii_in(i));
        
        dWni_Win=abs(Wni_l(i,:)-Win_l(i))./Win_l(i);
        Rw_ni_in_l(i)=wi_2(min(find(dWni_Win<=min(dWni_Win))));
        Fi_2_Rw_ni_in(i)=Fi_2(find(dWni_Win<=min(dWni_Win), 1 ));
        dWfn_Win=abs(Wfn_l(i,:)-Win_l(i))./Win_l(i);
        Rw_fn_in_l(i)=wf_1(find(dWfn_Win<=min(dWfn_Win), 1 ));
        Ff_1_Rw_fn_in(i)=Ff_1(find(dWfn_Win<=min(dWfn_Win), 1 ));
        dWif_Win=abs(Wif_l(i,:)-Win_l(i))./Win_l(i);
        Rw_if_in_l(i)=wf_2(find(dWif_Win<=min(dWif_Win), 1 ));
        Ff_2_Rw_if_in(i)=Ff_2(find(dWif_Win<=min(dWif_Win), 1 ));
        dWnf_Win=abs(Wnf_l(i,:)-Win_l(i))./Win_l(i);
        Rw_nf_in_l(i)=wf_2(find(dWnf_Win<=min(dWnf_Win), 1 ));
        Ff_2_Rw_nf_in(i)=Ff_2(find(dWnf_Win<=min(dWnf_Win), 1 ));
        Ff_2_minRw_if_in_nf_in(i)=Ff_2(min(find(dWif_Win<=min(dWif_Win), 1 ),min(find(dWnf_Win<=min(dWnf_Win)))));
        dWif_Wnf=abs(Wif_l(i,:)'-Wnf_l)./Wnf_l;
        Rw_if_nf_l(i)=wf_2(min(find(dWif_Wnf<=min(dWif_Wnf))));
        Ff_2_Rw_if_nf(i)=Ff_2(min(find(dWif_Wnf<=min(dWif_Wnf))));
        dWii_Win=abs(Wii_l(i,:)-Win_l(i))./Win_l(i);
        Rw_ii_in_l(i)=wi_2(min(find(dWii_Win<=min(dWii_Win))));
        Fi_2_Rw_ii_in(i)=Fi_2(min(find(dWii_Win<=min(dWii_Win))));
        
        
        derWin_l(i)=1./(r+di_1*(1-q_2)+di_1*q_2*Fi_2_Rw_ni_in(i)+lif_1*Ff_1_Rw_fn_in(i) ...
            +lnf_2*Ff_2_minRw_if_in_nf_in(i)+lni_2*Fi_2_Rw_ii_in(i));
    end
    
    for j=1:n
        % to construct Wnf_h
        dWin_Wnf=abs(Win_h-Wnf_h(j))./Wnf_h(j);
        Rw_in_nf_h(j)=wi_1(min(find(dWin_Wnf<=min(dWin_Wnf))));
        Fi_1_Rw_in_nf(j)=Fi_1(min(find(dWin_Wnf<=min(dWin_Wnf))));
        dWni_Wnf=abs(Wni_h-Wnf_h(j))./Wnf_h(j);
        Rw_ni_nf_h(j)=wi_2(min(find(dWni_Wnf<=min(dWni_Wnf))));
        Fi_2_Rw_ni_nf(j)=Fi_2(min(find(dWni_Wnf<=min(dWni_Wnf))));
        dWff_Wnf=abs(Wff_h(:,j)-Wnf_h(j))./Wnf_h(j);
        Rw_ff_nf_h(j)=wf_1(find(dWff_Wnf<=min(dWff_Wnf), 1 ));
        Ff_1_Rw_ff_nf(j)=Ff_1(find(dWff_Wnf<=min(dWff_Wnf), 1 ));
        dWfn_Wnf=abs(Wfn_h-Wnf_h(j))./Wnf_h(j);
        Rw_fn_nf_h(j)=wf_1(find(dWfn_Wnf<=min(dWfn_Wnf), 1 ));
        Ff_1_Rw_fn_nf(j)=Ff_1(min(find(dWfn_Wnf<=min(dWfn_Wnf))));
        Ff_1_minRw_ff_nf_fn_nf(j)=Ff_1(min(min(find(dWff_Wnf<=min(dWff_Wnf))),find(dWfn_Wnf<=min(dWfn_Wnf), 1 )));
        dWff_Wfn=abs(Wff_h(:,j)-Wfn_h)./Wfn_h;
        Rw_ff_fn_h(j)=wf_1(find(dWff_Wfn<=min(dWff_Wfn), 1 ));
        Ff_1_Rw_ff_fn(j)=Ff_1(find(dWff_Wfn<=min(dWff_Wfn), 1 ));
        dWif_Wnf=abs(Wif_h(:,j)-Wnf_h(j))./Wnf_h(j);
        Rw_if_nf_h(j)=wi_1(min(find(dWif_Wnf<=min(dWif_Wnf))));
        Fi_1_Rw_if_nf(j)=Fi_1(find(dWif_Wnf<=min(dWif_Wnf), 1 ));
        % to construct Wni_h
        dWin_Wni=abs(Win_h-Wni_h(j))./Wni_h(j);
        Rw_in_ni_h(j)=wi_1(find(dWin_Wni<=min(dWin_Wni), 1 ));
        Fi_1_Rw_in_ni(j)=Fi_1(find(dWin_Wni<=min(dWin_Wni), 1 ));
        dWnf_Wni=abs(Wnf_h-Wni_h(j))./Wni_h(j);
        Rw_nf_ni_h(j)=wf_2(min(find(dWnf_Wni<=min(dWnf_Wni))));
        Ff_2_Rw_nf_ni(j)=Ff_2(find(dWnf_Wni<=min(dWnf_Wni), 1 ));
        dWfi_Wni=abs(Wfi_h(:,j)-Wni_h(j))./Wni_h(j);
        Rw_fi_ni_h(j)=wf_1(min(find(dWfi_Wni<=min(dWfi_Wni))));
        Ff_1_Rw_fi_ni(j)=Ff_1(find(dWfi_Wni<=min(dWfi_Wni), 1 ));
        dWfn_Wni=abs(Wfn_h-Wni_h(j))./Wni_h(j);
        Rw_fn_ni_h(j)=wf_1(find(dWfn_Wni<=min(dWfn_Wni), 1 ));
        Ff_1_Rw_fn_ni(j)=Ff_1(find(dWfn_Wni<=min(dWfn_Wni), 1 ));
        Ff_1_minRw_fi_ni_fn_ni(j)=Ff_1(min(min(find(dWfi_Wni<=min(dWfi_Wni))),min(find(dWfn_Wni<=min(dWfn_Wni)))));
        dWfi_Wfn=abs(Wfi_h(:,j)-Wfn_h)./Wfn_h;
        Rw_fi_fn_h(j)=wf_1(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
        Ff_1_Rw_fi_fn(j)=Ff_1(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
        dWii_Wni=abs(Wii_h(:,j)-Wni_h(j))./Wni_h(j);
        Rw_ii_ni_h(j)=wi_1(min(find(dWii_Wni<=min(dWii_Wni))));
        Fi_1_Rw_ii_ni(j)=Fi_1(min(find(dWii_Wni<=min(dWii_Wni))));
        
        derWnf_h(j)=1./(r+df_2*(1-p_1)+df_2*p_1*Fi_1_Rw_in_nf(j)+lfi_2*Fi_2_Rw_ni_nf(j) ...
            +lnf_1*Ff_1_minRw_ff_nf_fn_nf(j)+lni_1*Fi_1_Rw_if_nf(j));
        
        derWni_h(j)=1./(r+di_2*(1-q_1)+di_2*q_1*Fi_1_Rw_in_ni(j)+lif_2*Ff_2_Rw_nf_ni(j) ...
            +lnf_1*Ff_1_minRw_fi_ni_fn_ni(j)+lni_1*Fi_1_Rw_ii_ni(j));
        
        dWin_Wnf=abs(Win_l-Wnf_l(j))./Wnf_l(j);
        Rw_in_nf_l(j)=wi_1(min(find(dWin_Wnf<=min(dWin_Wnf))));
        Fi_1_Rw_in_nf(j)=Fi_1(min(find(dWin_Wnf<=min(dWin_Wnf))));
        dWni_Wnf=abs(Wni_l-Wnf_l(j))./Wnf_l(j);
        Rw_ni_nf_l(j)=wi_2(min(find(dWni_Wnf<=min(dWni_Wnf))));
        Fi_2_Rw_ni_nf(j)=Fi_2(min(find(dWni_Wnf<=min(dWni_Wnf))));
        dWff_Wnf=abs(Wff_l(:,j)-Wnf_l(j))./Wnf_l(j);
        Rw_ff_nf_l(j)=wf_1(find(dWff_Wnf<=min(dWff_Wnf), 1 ));
        Ff_1_Rw_ff_nf(j)=Ff_1(find(dWff_Wnf<=min(dWff_Wnf), 1 ));
        dWfn_Wnf=abs(Wfn_l-Wnf_l(j))./Wnf_l(j);
        Rw_fn_nf_l(j)=wf_1(find(dWfn_Wnf<=min(dWfn_Wnf), 1 ));
        Ff_1_Rw_fn_nf(j)=Ff_1(min(find(dWfn_Wnf<=min(dWfn_Wnf))));
        Ff_1_minRw_ff_nf_fn_nf(j)=Ff_1(min(min(find(dWff_Wnf<=min(dWff_Wnf))),find(dWfn_Wnf<=min(dWfn_Wnf), 1 )));
        dWff_Wfn=abs(Wff_l(:,j)-Wfn_l)./Wfn_l;
        Rw_ff_fn_l(j)=wf_1(find(dWff_Wfn<=min(dWff_Wfn), 1 ));
        Ff_1_Rw_ff_fn(j)=Ff_1(find(dWff_Wfn<=min(dWff_Wfn), 1 ));
        dWif_Wnf=abs(Wif_l(:,j)-Wnf_l(j))./Wnf_l(j);
        Rw_if_nf_l(j)=wi_1(min(find(dWif_Wnf<=min(dWif_Wnf))));
        Fi_1_Rw_if_nf(j)=Fi_1(find(dWif_Wnf<=min(dWif_Wnf), 1 ));
        % to construct Wni_l
        dWin_Wni=abs(Win_l-Wni_l(j))./Wni_l(j);
        Rw_in_ni_l(j)=wi_1(find(dWin_Wni<=min(dWin_Wni), 1 ));
        Fi_1_Rw_in_ni(j)=Fi_1(find(dWin_Wni<=min(dWin_Wni), 1 ));
        dWnf_Wni=abs(Wnf_l-Wni_l(j))./Wni_l(j);
        Rw_nf_ni_l(j)=wf_2(min(find(dWnf_Wni<=min(dWnf_Wni))));
        Ff_2_Rw_nf_ni(j)=Ff_2(find(dWnf_Wni<=min(dWnf_Wni), 1 ));
        dWfi_Wni=abs(Wfi_l(:,j)-Wni_l(j))./Wni_l(j);
        Rw_fi_ni_l(j)=wf_1(min(find(dWfi_Wni<=min(dWfi_Wni))));
        Ff_1_Rw_fi_ni(j)=Ff_1(find(dWfi_Wni<=min(dWfi_Wni), 1 ));
        dWfn_Wni=abs(Wfn_l-Wni_l(j))./Wni_l(j);
        Rw_fn_ni_l(j)=wf_1(find(dWfn_Wni<=min(dWfn_Wni), 1 ));
        Ff_1_Rw_fn_ni(j)=Ff_1(find(dWfn_Wni<=min(dWfn_Wni), 1 ));
        Ff_1_minRw_fi_ni_fn_ni(j)=Ff_1(min(min(find(dWfi_Wni<=min(dWfi_Wni))),min(find(dWfn_Wni<=min(dWfn_Wni)))));
        dWfi_Wfn=abs(Wfi_l(:,j)-Wfn_l)./Wfn_l;
        Rw_fi_fn_l(j)=wf_1(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
        Ff_1_Rw_fi_fn(j)=Ff_1(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
        dWii_Wni=abs(Wii_l(:,j)-Wni_l(j))./Wni_l(j);
        Rw_ii_ni_l(j)=wi_1(min(find(dWii_Wni<=min(dWii_Wni))));
        Fi_1_Rw_ii_ni(j)=Fi_1(min(find(dWii_Wni<=min(dWii_Wni))));
        
        derWnf_l(j)=1./(r+df_2*(1-p_1)+df_2*p_1*Fi_1_Rw_in_nf(j)+lfi_2*Fi_2_Rw_ni_nf(j) ...
            +lnf_1*Ff_1_minRw_ff_nf_fn_nf(j)+lni_1*Fi_1_Rw_if_nf(j));
        
        derWni_l(j)=1./(r+di_2*(1-q_1)+di_2*q_1*Fi_1_Rw_in_ni(j)+lif_2*Ff_2_Rw_nf_ni(j) ...
            +lnf_1*Ff_1_minRw_fi_ni_fn_ni(j)+lni_1*Fi_1_Rw_ii_ni(j));
    end
    
    for i=1:n
        for j=1:n
            
            % to construct Wff_h
            dWif_Wff=abs(Wif_h(:,j)-Wff_h(i,j))./Wff_h(i,j);dWif_Wff=abs(Wif_h(:,n)-Wff_h(i,n))./Wff_h(i,n);
            Rw_if_ff_h(i,j)=wi_1(min(find(dWif_Wff<=min(dWif_Wff))));
            Fi_1_Rw_if_ff(i,j)=Fi_1(min(find(dWif_Wff<=min(dWif_Wff))));
            dWfi_Wff=abs(Wfi_h(i,:)-Wff_h(i,j))./Wff_h(i,j);dWfi_Wff=abs(Wfi_h(n,:)-Wff_h(n,j))./Wff_h(n,j);
            Rw_fi_ff_h(i,j)=wi_2(min(find(dWfi_Wff<=min(dWfi_Wff))));
            Fi_2_Rw_fi_ff(i,j)=Fi_2(min(find(dWfi_Wff<=min(dWfi_Wff))));
            
            derWff_h(i,j)=1./(r+df_1+df_2+lfi_1*Fi_1_Rw_if_ff(i,j) ...
                +lfi_2*Fi_2_Rw_fi_ff(i,j));
            
            dWif_Wff=abs(Wif_l(:,j)-Wff_l(i,j))./Wff_l(i,j);dWif_Wff=abs(Wif_l(:,n)-Wff_l(i,n))./Wff_l(i,n);
            Rw_if_ff_l(i,j)=wi_1(min(find(dWif_Wff<=min(dWif_Wff))));
            Fi_1_Rw_if_ff(i,j)=Fi_1(min(find(dWif_Wff<=min(dWif_Wff))));
            dWfi_Wff=abs(Wfi_l(i,:)-Wff_l(i,j))./Wff_l(i,j);dWfi_Wff=abs(Wfi_l(n,:)-Wff_l(n,j))./Wff_l(n,j);
            Rw_fi_ff_l(i,j)=wi_2(min(find(dWfi_Wff<=min(dWfi_Wff))));
            Fi_2_Rw_fi_ff(i,j)=Fi_2(min(find(dWfi_Wff<=min(dWfi_Wff))));
            
            derWff_l(i,j)=1./(r+df_1+df_2+lfi_1*Fi_1_Rw_if_ff(i,j) ...
                +lfi_2*Fi_2_Rw_fi_ff(i,j));
            
            % to construct Wii_h
            dWfi_Wii=abs(Wfi_h(:,j)-Wii_h(i,j))./Wii_h(i,j);dWfi_Wii=abs(Wfi_h(:,n)-Wii_h(i,n))./Wii_h(i,n);
            Rw_fi_ii_h(i,j)=wf_1(min(find(dWfi_Wii<=min(dWfi_Wii))));
            Ff_1_Rw_fi_ii(i,j)=Ff_1(find(dWfi_Wii<=min(dWfi_Wii), 1 ));
            dWfn_Wii=abs(Wfn_h-Wii_h(i,j))./Wii_h(i,j);
            Rw_fn_ii_h(i,j)=wf_1(min(find(dWfn_Wii<=min(dWfn_Wii))));
            Ff_1_Rw_fn_ii(i,j)=Ff_1(min(find(dWfn_Wii<=min(dWfn_Wii))));
            Ff_1_minRw_fi_ii_fn_ii(i,j)=Ff_1(min(min(find(dWfi_Wii<=min(dWfi_Wii))),find(dWfn_Wii<=min(dWfn_Wii), 1 )));
            dWfi_Wfn=abs(Wfi_h(:,j)-Wfn_h)./Wfn_h;
            Rw_fi_fn_h(j)=wf_1(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
            Ff_1_Rw_fi_fn(j)=Ff_1(find(dWfi_Wfn<=min(dWfi_Wfn), 1 ));
            dWif_Wii=abs(Wif_h(i,:)-Wii_h(i,j))./Wii_h(i,j);dWif_Wii=abs(Wif_h(n,:)-Wii_h(n,j))./Wii_h(n,j);
            Rw_if_ii_h(i,j)=wf_2(find(dWif_Wii<=min(dWif_Wii), 1 ));
            Ff_2_Rw_if_ii(i,j)=Ff_2(find(dWif_Wii<=min(dWif_Wii), 1 ));
            dWnf_Wii=abs(Wnf_h-Wii_h(i,j))./Wii_h(i,j);
            Rw_nf_ii_h(i,j)=wf_2(min(find(dWnf_Wii<=min(dWnf_Wii))));
            Ff_2_Rw_nf_ii(i,j)=Ff_2(find(dWnf_Wii<=min(dWnf_Wii), 1 ));
            Ff_2_minRw_if_ii_nf_ii(i,j)=Ff_2(min(find(dWif_Wii<=min(dWif_Wii), 1 ),min(find(dWnf_Wii<=min(dWnf_Wii)))));
            dWif_Wnf=abs(Wif_h(i,:)'-Wnf_h)./Wnf_h;
            Rw_if_nf_h(j)=wf_2(find(dWif_Wnf<=min(dWif_Wnf), 1 ));
            Ff_2_Rw_if_nf(j)=Ff_2(find(dWif_Wnf<=min(dWif_Wnf), 1 ));
            
            derWii_h(i,j)=1./(r+di_1+di_2+lif_1*Ff_1_minRw_fi_ii_fn_ii(i,j) ...
                +lif_2*Ff_2_minRw_if_ii_nf_ii(i,j));
            
            dWfi_Wii=abs(Wfi_l(:,j)-Wii_l(i,j))./Wii_l(i,j);dWfi_Wii=abs(Wfi_l(:,n)-Wii_l(i,n))./Wii_l(i,n);
            Rw_fi_ii_l(i,j)=wf_1(min(find(dWfi_Wii<=min(dWfi_Wii))));
            Ff_1_Rw_fi_ii(i,j)=Ff_1(find(dWfi_Wii<=min(dWfi_Wii), 1 ));
            dWfn_Wii=abs(Wfn_l-Wii_l(i,j))./Wii_l(i,j);
            Rw_fn_ii_l(i,j)=wf_1(min(find(dWfn_Wii<=min(dWfn_Wii))));
            Ff_1_Rw_fn_ii(i,j)=Ff_1(min(find(dWfn_Wii<=min(dWfn_Wii))));
            Ff_1_minRw_fi_ii_fn_ii(i,j)=Ff_1(min(min(find(dWfi_Wii<=min(dWfi_Wii))),find(dWfn_Wii<=min(dWfn_Wii), 1 )));
            dWfi_Wfn=abs(Wfi_l(:,j)-Wfn_l)./Wfn_l;
            Rw_fi_fn_l(j)=wf_1(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
            Ff_1_Rw_fi_fn(j)=Ff_1(find(dWfi_Wfn<=min(dWfi_Wfn), 1 ));
            dWif_Wii=abs(Wif_l(i,:)-Wii_l(i,j))./Wii_l(i,j);dWif_Wii=abs(Wif_l(n,:)-Wii_l(n,j))./Wii_l(n,j);
            Rw_if_ii_l(i,j)=wf_2(find(dWif_Wii<=min(dWif_Wii), 1 ));
            Ff_2_Rw_if_ii(i,j)=Ff_2(find(dWif_Wii<=min(dWif_Wii), 1 ));
            dWnf_Wii=abs(Wnf_l-Wii_l(i,j))./Wii_l(i,j);
            Rw_nf_ii_l(i,j)=wf_2(min(find(dWnf_Wii<=min(dWnf_Wii))));
            Ff_2_Rw_nf_ii(i,j)=Ff_2(find(dWnf_Wii<=min(dWnf_Wii), 1 ));
            Ff_2_minRw_if_ii_nf_ii(i,j)=Ff_2(min(find(dWif_Wii<=min(dWif_Wii), 1 ),min(find(dWnf_Wii<=min(dWnf_Wii)))));
            dWif_Wnf=abs(Wif_l(i,:)'-Wnf_l)./Wnf_l;
            Rw_if_nf_l(j)=wf_2(find(dWif_Wnf<=min(dWif_Wnf), 1 ));
            Ff_2_Rw_if_nf(j)=Ff_2(find(dWif_Wnf<=min(dWif_Wnf), 1 ));
            
            derWii_l(i,j)=1./(r+di_1+di_2+lif_1*Ff_1_minRw_fi_ii_fn_ii(i,j) ...
                +lif_2*Ff_2_minRw_if_ii_nf_ii(i,j));
            
            % to construct Wfi_h
            dWii_Wfi=abs(Wii_h(:,j)-Wfi_h(i,j))./Wfi_h(i,j);dWii_Wfi=abs(Wii_h(:,n)-Wfi_h(i,n))./Wfi_h(i,n);
            Rw_ii_fi_h(i,j)=wi_1(min(find(dWii_Wfi<=min(dWii_Wfi))));
            Fi_1_Rw_ii_fi(i,j)=Fi_1(find(dWii_Wfi<=min(dWii_Wfi), 1 ));
            dWff_Wfi=abs(Wff_h(i,:)-Wfi_h(i,j))./Wfi_h(i,j);dWff_Wfi=abs(Wff_h(n,:)-Wfi_h(n,j))./Wfi_h(n,j);
            Rw_ff_fi_h(i,j)=wf_2(min(find(dWff_Wfi<=min(dWff_Wfi))));
            Ff_2_Rw_ff_fi(i,j)=Ff_2(min(find(dWff_Wfi<=min(dWff_Wfi))));
            dWnf_Wfi=abs(Wnf_h-Wfi_h(i,j))./Wfi_h(i,j);
            Rw_nf_fi_h(i,j)=wf_2(find(dWnf_Wfi<=min(dWnf_Wfi), 1 ));
            Ff_2_Rw_nf_fi(i,j)=Ff_2(min(find(dWnf_Wfi<=min(dWnf_Wfi))));
            Ff_2_minRw_ff_fi_nf_fi(i,j)=Ff_2(min(min(find(dWff_Wfi<=min(dWff_Wfi))),min(find(dWnf_Wfi<=min(dWnf_Wfi)))));
            dWff_Wnf=abs(Wff_h(i,:)'-Wnf_h)./Wnf_h;
            Rw_ff_nf_h(i)=wf_2(find(dWff_Wnf<=min(dWff_Wnf), 1 ));
            Ff_2_Rw_ff_nf(i)=Ff_2(find(dWff_Wnf<=min(dWff_Wnf), 1 ));
            
            derWfi_h(i,j)=1./(r+df_1+di_2+lfi_1*Fi_1_Rw_ii_fi(i,j) ...
                +lif_2*Ff_2_minRw_ff_fi_nf_fi(i,j));
            
            dWii_Wfi=abs(Wii_l(:,j)-Wfi_l(i,j))./Wfi_l(i,j);dWii_Wfi=abs(Wii_l(:,n)-Wfi_l(i,n))./Wfi_l(i,n);
            Rw_ii_fi_l(i,j)=wi_1(min(find(dWii_Wfi<=min(dWii_Wfi))));
            Fi_1_Rw_ii_fi(i,j)=Fi_1(find(dWii_Wfi<=min(dWii_Wfi), 1 ));
            dWff_Wfi=abs(Wff_l(i,:)-Wfi_l(i,j))./Wfi_l(i,j);dWff_Wfi=abs(Wff_l(n,:)-Wfi_l(n,j))./Wfi_l(n,j);
            Rw_ff_fi_l(i,j)=wf_2(min(find(dWff_Wfi<=min(dWff_Wfi))));
            Ff_2_Rw_ff_fi(i,j)=Ff_2(min(find(dWff_Wfi<=min(dWff_Wfi))));
            dWnf_Wfi=abs(Wnf_l-Wfi_l(i,j))./Wfi_l(i,j);
            Rw_nf_fi_l(i,j)=wf_2(find(dWnf_Wfi<=min(dWnf_Wfi), 1 ));
            Ff_2_Rw_nf_fi(i,j)=Ff_2(min(find(dWnf_Wfi<=min(dWnf_Wfi))));
            Ff_2_minRw_ff_fi_nf_fi(i,j)=Ff_2(min(min(find(dWff_Wfi<=min(dWff_Wfi))),min(find(dWnf_Wfi<=min(dWnf_Wfi)))));
            dWff_Wnf=abs(Wff_l(i,:)'-Wnf_l)./Wnf_l;
            Rw_ff_nf_l(i)=wf_2(find(dWff_Wnf<=min(dWff_Wnf), 1 ));
            Ff_2_Rw_ff_nf(i)=Ff_2(find(dWff_Wnf<=min(dWff_Wnf), 1 ));
            
            derWfi_l(i,j)=1./(r+df_1+di_2+lfi_1*Fi_1_Rw_ii_fi(i,j) ...
                +lif_2*Ff_2_minRw_ff_fi_nf_fi(i,j));
            
            % to construct Wif_h
            dWff_Wif=abs(Wff_h(:,j)-Wif_h(i,j))./Wif_h(i,j);dWff_Wif=abs(Wff_h(:,n)-Wif_h(i,n))./Wif_h(i,n);
            Rw_ff_if_h(i,j)=wf_1(min(find(dWff_Wif<=min(dWff_Wif))));
            Ff_1_Rw_ff_if(i,j)=Ff_1(min(find(dWff_Wif<=min(dWff_Wif))));
            dWfn_Wif=abs(Wfn_h-Wif_h(i,j))./Wif_h(i,j);
            Rw_fn_if_h(i,j)=wf_1(min(find(dWfn_Wif<=min(dWfn_Wif))));
            Ff_1_Rw_fn_if(i,j)=Ff_1(find(dWfn_Wif<=min(dWfn_Wif), 1 ));
            Ff_1_minRw_ff_if_fn_if(i,j)=Ff_1(min(find(dWff_Wif<=min(dWff_Wif), 1 ),find(dWfn_Wif<=min(dWfn_Wif), 1 )));
            dWff_Wfn=abs(Wff_h(:,j)-Wfn_h)./Wfn_h;
            Rw_ff_fn_h(j)=wf_1(find(dWff_Wfn<=min(dWff_Wfn), 1 ));
            Ff_1_Rw_ff_fn(j)=Ff_1(min(find(dWff_Wfn<=min(dWff_Wfn))));
            dWii_Wif=abs(Wii_h(i,:)-Wif_h(i,j))./Wif_h(i,j);dWii_Wif=abs(Wii_h(n,:)-Wif_h(n,j))./Wif_h(n,j);
            Rw_ii_if_h(i,j)=wi_2(min(find(dWii_Wif<=min(dWii_Wif))));
            Fi_2_Rw_ii_if(i,j)=Fi_2(find(dWii_Wif<=min(dWii_Wif), 1 ));
            
            derWif_h(i,j)=1./(r+di_1+df_2+lif_1*Ff_1_minRw_ff_if_fn_if(i,j) ...
                +lfi_2*Fi_2_Rw_ii_if(i,j));
            dWff_Wif=abs(Wff_l(:,j)-Wif_l(i,j))./Wif_l(i,j);dWff_Wif=abs(Wff_l(:,n)-Wif_l(i,n))./Wif_l(i,n);
            Rw_ff_if_l(i,j)=wf_1(min(find(dWff_Wif<=min(dWff_Wif))));
            Ff_1_Rw_ff_if(i,j)=Ff_1(min(find(dWff_Wif<=min(dWff_Wif))));
            dWfn_Wif=abs(Wfn_l-Wif_l(i,j))./Wif_l(i,j);
            Rw_fn_if_l(i,j)=wf_1(min(find(dWfn_Wif<=min(dWfn_Wif))));
            Ff_1_Rw_fn_if(i,j)=Ff_1(find(dWfn_Wif<=min(dWfn_Wif), 1 ));
            Ff_1_minRw_ff_if_fn_if(i,j)=Ff_1(min(find(dWff_Wif<=min(dWff_Wif), 1 ),find(dWfn_Wif<=min(dWfn_Wif), 1 )));
            dWff_Wfn=abs(Wff_l(:,j)-Wfn_l)./Wfn_l;
            Rw_ff_fn_l(j)=wf_1(find(dWff_Wfn<=min(dWff_Wfn), 1 ));
            Ff_1_Rw_ff_fn(j)=Ff_1(min(find(dWff_Wfn<=min(dWff_Wfn))));
            dWii_Wif=abs(Wii_l(i,:)-Wif_l(i,j))./Wif_l(i,j);dWii_Wif=abs(Wii_l(n,:)-Wif_l(n,j))./Wif_l(n,j);
            Rw_ii_if_l(i,j)=wi_2(min(find(dWii_Wif<=min(dWii_Wif))));
            Fi_2_Rw_ii_if(i,j)=Fi_2(find(dWii_Wif<=min(dWii_Wif), 1 ));
            
            derWif_l(i,j)=1./(r+di_1+df_2+lif_1*Ff_1_minRw_ff_if_fn_if(i,j) ...
                +lfi_2*Fi_2_Rw_ii_if(i,j));
            
            %here
            % to construct Wnn
            dWfn_Wnn=abs(Wfn_l-Wnn_l)./Wnn_l;
            Rw_fn_nn_l=wf_1(find(dWfn_Wnn<=min(dWfn_Wnn), 1 ));
            %Ff_1_Rw_fn_nn=Ff_1(min(find(dWfn_Wnn<=min(dWfn_Wnn))));
            
            dWnf_Wnn=abs(Wnf_l-Wnn_l)./Wnn_l;
            Rw_nf_nn_l=wf_2(min(find(dWnf_Wnn<=min(dWnf_Wnn))));
            %Ff_2_Rw_nf_nn=Ff_2(min(find(dWnf_Wnn<=min(dWnf_Wnn))));
            
            dWin_Wnn=abs(Win_l-Wnn_l)./Wnn_l;
            Rw_in_nn_l=wi_1(find(dWin_Wnn<=min(dWin_Wnn), 1 ));
            %Fi_1_Rw_in_nn=Fi_1(min(find(dWin_Wnn<=min(dWin_Wnn))));
            
            dWni_Wnn=abs(Wni_l-Wnn_l)./Wnn_l;
            Rw_ni_nn_l=wi_2(find(dWni_Wnn<=min(dWni_Wnn), 1 ));
            %Fi_2_Rw_ni_nn=Fi_2(find(dWni_Wnn<=min(dWni_Wnn), 1 ));
            
            dWfn_Wnn=abs(Wfn_h-Wnn_h)./Wnn_h;
            Rw_fn_nn_h=wf_1(find(dWfn_Wnn<=min(dWfn_Wnn), 1 ));
            %Ff_1_Rw_fn_nn=Ff_1(min(find(dWfn_Wnn<=min(dWfn_Wnn))));
            
            dWnf_Wnn=abs(Wnf_h-Wnn_h)./Wnn_h;
            Rw_nf_nn_h=wf_2(min(find(dWnf_Wnn<=min(dWnf_Wnn))));
            %Ff_2_Rw_nf_nn=Ff_2(min(find(dWnf_Wnn<=min(dWnf_Wnn))));
            
            dWin_Wnn=abs(Win_h-Wnn_h)./Wnn_h;
            Rw_in_nn_h=wi_1(find(dWin_Wnn<=min(dWin_Wnn), 1 ));
            %Fi_1_Rw_in_nn=Fi_1(min(find(dWin_Wnn<=min(dWin_Wnn))));
            
            dWni_Wnn=abs(Wni_h-Wnn_h)./Wnn_h;
            Rw_ni_nn_h=wi_2(find(dWni_Wnn<=min(dWni_Wnn), 1 ));
            %Fi_2_Rw_ni_nn=Fi_2(find(dWni_Wnn<=min(dWni_Wnn), 1 ));
            
            
        end
    end
    
    %here
    % solve integral by parts
    
    % Wfn_h
    for o=1:n
        dWff_Wfn=abs(Wff_h(o,:)-Wfn_h(o))./Wfn_h(o);
        Rw_ff_fn_h(o)=wf_2(min(find(dWff_Wfn<=min(dWff_Wfn))));
        dWff_Wnf=abs(Wff_h(o,:)'-Wnf_h)./Wnf_h;
        Rw_ff_nf_h(o)=wf_2(min(find(dWff_Wnf<=min(dWff_Wnf))));
        dWfi_Wfn=abs(Wfi_h(o,:)-Wfn_h(o))./Wfn_h(o);
        Rw_fi_fn_h(o)=wi_2(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
        
        dWff_Wfn=abs(Wff_l(o,:)-Wfn_l(o))./Wfn_l(o);
        Rw_ff_fn_l(o)=wf_2(min(find(dWff_Wfn<=min(dWff_Wfn))));
        dWff_Wnf=abs(Wff_l(o,:)'-Wnf_l)./Wnf_l;
        Rw_ff_nf_l(o)=wf_2(min(find(dWff_Wnf<=min(dWff_Wnf))));
        dWfi_Wfn=abs(Wfi_l(o,:)-Wfn_l(o))./Wfn_l(o);
        Rw_fi_fn_l(o)=wi_2(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
    end
    for i=1:n
        
        intmaxWni_Wfn_h(i)=sum((wi_2>Rw_ni_fn_h(i)).*Fi_2.*derWni_h'.*wht_i_2);
        intmaxWin_Wfn_h(i)=sum((wi_1>Rw_in_fn_h(i)).*Fi_1.*derWin_h'.*wht_i_1);
        intmaxWfi_Wfn_h(i)=sum((wi_2>Rw_fi_fn_h(i)).*Fi_2.*derWfi_h(i,:)'.*wht_i_2);
        %intmaxWff_Wnf_Wfn(i)=sum((wf_2>Rw_ff_fn(i)).*Ff_2.*derWff(i,:)'.*wht_f_2);
        intmaxWff_Wnf_Wfn_h(i)=max(0,sum((wf_2<Rw_ff_nf_h(i)).*Ff_2.*derWff_h(i,:)'.*wht_f_2+(wf_2>=Rw_ff_nf_h(i)).*Ff_2.*derWnf_h'.*wht_f_2));
        
        intmaxWni_Wfn_l(i)=sum((wi_2>Rw_ni_fn_l(i)).*Fi_2.*derWni_l'.*wht_i_2);
        intmaxWin_Wfn_l(i)=sum((wi_1>Rw_in_fn_l(i)).*Fi_1.*derWin_l'.*wht_i_1);
        intmaxWfi_Wfn_l(i)=sum((wi_2>Rw_fi_fn_l(i)).*Fi_2.*derWfi_l(i,:)'.*wht_i_2);
        %intmaxWff_Wnf_Wfn(i)=sum((wf_2>Rw_ff_fn(i)).*Ff_2.*derWff(i,:)'.*wht_f_2);
        intmaxWff_Wnf_Wfn_l(i)=max(0,sum((wf_2<Rw_ff_nf_l(i)).*Ff_2.*derWff_l(i,:)'.*wht_f_2+(wf_2>=Rw_ff_nf_l(i)).*Ff_2.*derWnf_l'.*wht_f_2));
        
    end
    % Win_h
    
    for o=1:n
        dWif_Wnf=abs(Wif_h(o,:)'-Wnf_h)./Wnf_h;
        Rw_if_nf_h(o)=wf_2(min(find(dWif_Wnf<=min(dWif_Wnf))));
        
        dWif_Wnf=abs(Wif_l(o,:)'-Wnf_l)./Wnf_l;
        Rw_if_nf_l(o)=wf_2(min(find(dWif_Wnf<=min(dWif_Wnf))));
    end
    
    for i=1:n
        
        intmaxWni_Win_h(i)=sum((wi_2>Rw_ni_in_h(i)).*Fi_2.*derWni_h'.*wht_i_2);
        intmaxWfn_Win_h(i)=sum((wf_1>Rw_fn_in_h(i)).*Ff_1.*derWfn_h'.*wht_f_1);
        intmaxWii_Win_h(i)=sum((wi_2>Rw_ii_in_h(i)).*Fi_2.*derWii_h(i,:)'.*wht_i_2);
        %intmaxWif_Wnf_Win(i)=sum((wf_2>Rw_if_in(i)).*Ff_2.*derWif(i,:)'.*wht_f_2);
        intmaxWif_Wnf_Win_h(i)=max(0,sum((wf_2<Rw_if_nf_h(i)).*Ff_2.*derWif_h(i,:)'.*wht_f_2+(wf_2>=Rw_if_nf_h(i)).*Ff_2.*derWnf_h'.*wht_f_2));
        
        intmaxWni_Win_l(i)=sum((wi_2>Rw_ni_in_l(i)).*Fi_2.*derWni_l'.*wht_i_2);
        intmaxWfn_Win_l(i)=sum((wf_1>Rw_fn_in_l(i)).*Ff_1.*derWfn_l'.*wht_f_1);
        intmaxWii_Win_l(i)=sum((wi_2>Rw_ii_in_l(i)).*Fi_2.*derWii_l(i,:)'.*wht_i_2);
        %intmaxWif_Wnf_Win(i)=sum((wf_2>Rw_if_in(i)).*Ff_2.*derWif(i,:)'.*wht_f_2);
        intmaxWif_Wnf_Win_l(i)=max(0,sum((wf_2<Rw_if_nf_l(i)).*Ff_2.*derWif_l(i,:)'.*wht_f_2+(wf_2>=Rw_if_nf_l(i)).*Ff_2.*derWnf_l'.*wht_f_2));
        
    end
    
    %here
    % Wnf_h
    for o=1:n
        dWff_Wfn=abs(Wff_h(:,o)-Wfn_h)./Wfn_h;
        Rw_ff_fn_h(o)=wf_1(min(find(dWff_Wfn<=min(dWff_Wfn))));
        dWff_Wnf=abs(Wff_h(:,o)-Wnf_h(o))./Wnf_h(o);
        Rw_ff_nf_h(o)=wf_1(min(find(dWff_Wnf<=min(dWff_Wnf))));
        dWif_Wnf=abs(Wif_h(:,o)-Wnf_h(o))./Wnf_h(o);
        Rw_if_nf_h(o)=wi_1(min(find(dWif_Wnf<=min(dWif_Wnf))));
        
        dWff_Wfn=abs(Wff_l(:,o)-Wfn_l)./Wfn_l;
        Rw_ff_fn_l(o)=wf_1(min(find(dWff_Wfn<=min(dWff_Wfn))));
        dWff_Wnf=abs(Wff_l(:,o)-Wnf_l(o))./Wnf_l(o);
        Rw_ff_nf_l(o)=wf_1(min(find(dWff_Wnf<=min(dWff_Wnf))));
        dWif_Wnf=abs(Wif_l(:,o)-Wnf_l(o))./Wnf_l(o);
        Rw_if_nf_l(o)=wi_1(min(find(dWif_Wnf<=min(dWif_Wnf))));
    end
    
    for i=1:n
        
        intmaxWin_Wnf_h(i)=sum((wi_1>Rw_in_nf_h(i)).*Fi_1.*derWin_h'.*wht_i_1);
        intmaxWni_Wnf_h(i)=sum((wi_2>Rw_ni_nf_h(i)).*Fi_2.*derWni_h'.*wht_i_2);
        intmaxWif_Wnf_h(i)=sum((wi_1>Rw_if_nf_h(i)).*Fi_1.*derWif_h(:,i).*wht_i_1);
        %intmaxWff_Wfn_Wnf(i)=sum((wf_1>Rw_ff_nf(i)).*Ff_1.*derWff(:,i).*wht_f_1);
        intmaxWff_Wfn_Wnf_h(i)=max(0,sum((wf_1<Rw_ff_fn_h(i)).*Ff_1.*derWff_h(:,i).*wht_f_1+(wf_1>=Rw_ff_fn_h(i)).*Ff_1.*derWfn_h'.*wht_f_1));
        
        intmaxWin_Wnf_l(i)=sum((wi_1>Rw_in_nf_l(i)).*Fi_1.*derWin_l'.*wht_i_1);
        intmaxWni_Wnf_l(i)=sum((wi_2>Rw_ni_nf_l(i)).*Fi_2.*derWni_l'.*wht_i_2);
        intmaxWif_Wnf_l(i)=sum((wi_1>Rw_if_nf_l(i)).*Fi_1.*derWif_l(:,i).*wht_i_1);
        %intmaxWff_Wfn_Wnf(i)=sum((wf_1>Rw_ff_nf(i)).*Ff_1.*derWff(:,i).*wht_f_1);
        intmaxWff_Wfn_Wnf_l(i)=max(0,sum((wf_1<Rw_ff_fn_l(i)).*Ff_1.*derWff_l(:,i).*wht_f_1+(wf_1>=Rw_ff_fn_l(i)).*Ff_1.*derWfn_l'.*wht_f_1));
        
    end
    
    % Wni_h
    for o=1:n
        dWfi_Wfn=abs(Wfi_h(:,o)-Wfn_h)./Wfn_h;
        Rw_fi_fn_h(o)=wf_1(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
        
        dWfi_Wfn=abs(Wfi_l(:,o)-Wfn_l)./Wfn_l;
        Rw_fi_fn_l(o)=wf_1(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
    end
    
    for i=1:n
        
        intmaxWin_Wni_h(i)=sum((wi_1>Rw_in_ni_h(i)).*Fi_1.*derWin_h'.*wht_i_1);
        intmaxWnf_Wni_h(i)=sum((wf_2>Rw_nf_ni_h(i)).*Ff_2.*derWnf_h'.*wht_f_2);
        intmaxWii_Wni_h(i)=sum((wi_1>Rw_ii_ni_h(i)).*Fi_1.*derWii_h(:,i).*wht_i_1);
        %intmaxWfi_Wfn_Wni(i)=sum((wf_1>Rw_fi_ni(i)).*Ff_1.*derWfi(:,i).*wht_f_1);
        intmaxWfi_Wfn_Wni_h(i)=max(0,sum((wf_1<Rw_fi_fn_h(i)).*Ff_1.*derWfi_h(:,i).*wht_f_1+(wf_1>=Rw_fi_fn_h(i)).*Ff_1.*derWfn_h'.*wht_f_1));
        
        intmaxWin_Wni_l(i)=sum((wi_1>Rw_in_ni_l(i)).*Fi_1.*derWin_l'.*wht_i_1);
        intmaxWnf_Wni_l(i)=sum((wf_2>Rw_nf_ni_l(i)).*Ff_2.*derWnf_l'.*wht_f_2);
        intmaxWii_Wni_l(i)=sum((wi_1>Rw_ii_ni_l(i)).*Fi_1.*derWii_l(:,i).*wht_i_1);
        %intmaxWfi_Wfn_Wni(i)=sum((wf_1>Rw_fi_ni(i)).*Ff_1.*derWfi(:,i).*wht_f_1);
        intmaxWfi_Wfn_Wni_l(i)=max(0,sum((wf_1<Rw_fi_fn_l(i)).*Ff_1.*derWfi_l(:,i).*wht_f_1+(wf_1>=Rw_fi_fn_l(i)).*Ff_1.*derWfn_l'.*wht_f_1));
    end
    
    
    for i=1:n
        for j=1:n
            
            % Wff_h
            
            intmaxWif_Wff_h(i,j)=sum((wi_1>Rw_if_ff_h(i,j)).*Fi_1.*derWif_h(:,j).*wht_i_1);
            intmaxWfi_Wff_h(i,j)=sum((wi_2>Rw_fi_ff_h(i,j)).*Fi_2.*derWfi_h(i,:)'.*wht_i_2);
            
            intmaxWif_Wff_l(i,j)=sum((wi_1>Rw_if_ff_l(i,j)).*Fi_1.*derWif_l(:,j).*wht_i_1);
            intmaxWfi_Wff_l(i,j)=sum((wi_2>Rw_fi_ff_l(i,j)).*Fi_2.*derWfi_l(i,:)'.*wht_i_2);
        end
    end
    
    % Wii_h
    for o=1:n
        dWfi_Wfn=abs(Wfi_h(:,o)-Wfn_h)./Wfn_h;
        Rw_fi_fn_h(o)=wf_1(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
        dWif_Wnf=abs(Wif_h(o,:)'-Wnf_h)./Wnf_h;
        Rw_if_nf_h(o)=wf_2(min(find(dWif_Wnf<=min(dWif_Wnf))));
        
        dWfi_Wfn=abs(Wfi_l(:,o)-Wfn_l)./Wfn_l;
        Rw_fi_fn_l(o)=wf_1(min(find(dWfi_Wfn<=min(dWfi_Wfn))));
        dWif_Wnf=abs(Wif_l(o,:)'-Wnf_l)./Wnf_l;
        Rw_if_nf_l(o)=wf_2(min(find(dWif_Wnf<=min(dWif_Wnf))));
    end
    
    for i=1:n
        for j=1:n
            
            % intmaxWfi_Wfn_Wii(i,j)=sum((wf_1>Rw_fi_ii(i,j)).*Ff_1.*derWfi(:,j).*wht_f_1);
            % intmaxWif_Wnf_Wii(i,j)=sum((wf_2>Rw_if_ii(i,j)).*Ff_2.*derWif(i,:)'.*wht_f_2);
            
            intmaxWfi_Wfn_Wii_h(i,j)=max(0,sum((wf_1<Rw_fi_fn_h(j)).*Ff_1.*derWfi_h(:,j).*wht_f_1+(wf_1>=Rw_fi_fn_h(j)).*Ff_1.*derWfn_h'.*wht_f_1));
            intmaxWif_Wnf_Wii_h(i,j)=max(0,sum((wf_2<Rw_if_nf_h(i)).*Ff_2.*derWif_h(i,:)'.*wht_f_2+(wf_2>=Rw_if_nf_h(i)).*Ff_2.*derWnf_h'.*wht_f_2));
            
            intmaxWfi_Wfn_Wii_l(i,j)=max(0,sum((wf_1<Rw_fi_fn_l(j)).*Ff_1.*derWfi_l(:,j).*wht_f_1+(wf_1>=Rw_fi_fn_l(j)).*Ff_1.*derWfn_l'.*wht_f_1));
            intmaxWif_Wnf_Wii_l(i,j)=max(0,sum((wf_2<Rw_if_nf_l(i)).*Ff_2.*derWif_l(i,:)'.*wht_f_2+(wf_2>=Rw_if_nf_l(i)).*Ff_2.*derWnf_l'.*wht_f_2));
            
        end
    end
    
    
    % Wfi_h
    for o=1:n
        dWff_Wnf=abs(Wff_h(o,:)'-Wnf_h)./Wnf_h;
        Rw_ff_nf_h(o)=wf_2(min(find(dWff_Wnf<=min(dWff_Wnf))));
        
        dWff_Wnf=abs(Wff_l(o,:)'-Wnf_l)./Wnf_l;
        Rw_ff_nf_l(o)=wf_2(min(find(dWff_Wnf<=min(dWff_Wnf))));
    end
    
    
    for i=1:n
        for j=1:n
            intmaxWii_Wfi_h(i,j)=sum((wi_1>Rw_ii_fi_h(i,j)).*Fi_1.*derWii_h(:,j).*wht_i_1);
            % intmaxWff_Wnf_Wfi(i,j)=sum((wf_2>Rw_ff_fi(i,j)).*Ff_2.*derWff(i,:)'.*wht_f_2);
            intmaxWff_Wnf_Wfi_h(i,j)=max(0,sum((wf_2<Rw_ff_nf_h(i)).*Ff_2.*derWff_h(i,:)'.*wht_f_2+(wf_2>=Rw_ff_nf_h(i)).*Ff_2.*derWnf_h'.*wht_f_2));
            
            intmaxWii_Wfi_l(i,j)=sum((wi_1>Rw_ii_fi_l(i,j)).*Fi_1.*derWii_l(:,j).*wht_i_1);
            % intmaxWff_Wnf_Wfi(i,j)=sum((wf_2>Rw_ff_fi(i,j)).*Ff_2.*derWff(i,:)'.*wht_f_2);
            intmaxWff_Wnf_Wfi_l(i,j)=max(0,sum((wf_2<Rw_ff_nf_l(i)).*Ff_2.*derWff_l(i,:)'.*wht_f_2+(wf_2>=Rw_ff_nf_l(i)).*Ff_2.*derWnf_l'.*wht_f_2));
        end
    end
    
    % Wif_h
    for o=1:n
        dWff_Wfn_h=abs(Wff_h(:,o)-Wfn_h)./Wfn_h;
        Rw_ff_fn_h(o)=wf_1(min(find(dWff_Wfn_h<=min(dWff_Wfn_h))));
        
        dWff_Wfn_l=abs(Wff_l(:,o)-Wfn_l)./Wfn_l;
        Rw_ff_fn_l(o)=wf_1(min(find(dWff_Wfn_l<=min(dWff_Wfn_l))));
    end
    
    for i=1:n
        for j=1:n
            
            intmaxWii_Wif_h(i,j)=sum((wi_2>Rw_ii_if_h(i,j)).*Fi_2.*derWii_h(i,:)'.*wht_i_2);
            % intmaxWff_Wfn_Wif(i,j)=sum((wf_1>Rw_ff_if(i,j)).*Ff_1.*derWff(:,j).*wht_f_1);
            intmaxWff_Wfn_Wif_h(i,j)=max(0,sum((wf_1<Rw_ff_fn_h(j)).*Ff_1.*derWff_h(:,j).*wht_f_1+(wf_1>=Rw_ff_fn_h(j)).*Ff_1.*derWfn_h'.*wht_f_1));
            
            intmaxWii_Wif_l(i,j)=sum((wi_2>Rw_ii_if_l(i,j)).*Fi_2.*derWii_l(i,:)'.*wht_i_2);
            % intmaxWff_Wfn_Wif(i,j)=sum((wf_1>Rw_ff_if(i,j)).*Ff_1.*derWff(:,j).*wht_f_1);
            intmaxWff_Wfn_Wif_l(i,j)=max(0,sum((wf_1<Rw_ff_fn_l(j)).*Ff_1.*derWff_l(:,j).*wht_f_1+(wf_1>=Rw_ff_fn_l(j)).*Ff_1.*derWfn_l'.*wht_f_1));
        end
    end
    
    % Wnn_h
    
    intmaxWfn_Wnn_h=sum((wf_1>Rw_fn_nn_h).*Ff_1.*derWfn_h'.*wht_f_1);
    intmaxWin_Wnn_h=sum((wi_1>Rw_in_nn_h).*Fi_1.*derWin_h'.*wht_i_1);
    intmaxWnf_Wnn_h=sum((wf_2>Rw_nf_nn_h).*Ff_2.*derWnf_h'.*wht_f_2);
    intmaxWni_Wnn_h=sum((wi_2>Rw_ni_nn_h).*Fi_2.*derWni_h'.*wht_i_2);
    
    intmaxWfn_Wnn_l=sum((wf_1>Rw_fn_nn_l).*Ff_1.*derWfn_l'.*wht_f_1);
    intmaxWin_Wnn_l=sum((wi_1>Rw_in_nn_l).*Fi_1.*derWin_l'.*wht_i_1);
    intmaxWnf_Wnn_l=sum((wf_2>Rw_nf_nn_l).*Ff_2.*derWnf_l'.*wht_f_2);
    intmaxWni_Wnn_l=sum((wi_2>Rw_ni_nn_l).*Fi_2.*derWni_l'.*wht_i_2);
    
    
    
    %%%%% Value functions
    % solve in terms of F, reservation wages and lambdas
    
    % Wfn_h(w1) and Win_h(w1)
    for i=1:n
               
        eWfn_h(i)=(1/(r+df_1*(1-p_2)))*((1-exp(-theta*(wf_1(i)+b2_h)))/theta+a_h+df_1*(1-p_2)*Wnn_h ...
            +df_1*p_2*intmaxWni_Wfn_h(i) ...
            +lfi_1*intmaxWin_Wfn_h(i) ...
            +lnf_2*intmaxWff_Wnf_Wfn_h(i) ...
            +lni_2*intmaxWfi_Wfn_h(i))+...
            nu_hl*(1/(r+df_1*(1-p_2)))*(Wfn_l(i)-Wfn_h(i))+(1-nu_hl)*(1/(r+df_1*(1-p_2)))*0-Wfn_h(i);
       
                
        eWfn_l(i)=(1/(r+df_1*(1-p_2)))*((1-exp(-theta*(wf_1(i)+b2_l)))/theta+a_l+df_1*(1-p_2)*Wnn_l ...
            +df_1*p_2*intmaxWni_Wfn_l(i) ...
            +lfi_1*intmaxWin_Wfn_l(i) ...
            +lnf_2*intmaxWff_Wnf_Wfn_l(i) ...
            +lni_2*intmaxWfi_Wfn_l(i))+...
            nu_lh*(1/(r+df_1*(1-p_2)))*(Wfn_h(i)-Wfn_l(i))+(1-nu_lh)*(1/(r+df_1*(1-p_2)))*0-Wfn_l(i);
 
        eWin_h(i)=(1/(r+di_1*(1-q_2)))*((1-exp(-theta*(wi_1(i)+b2_h)))/theta+gamma+di_1*(1-q_2)*Wnn_h ...
            +di_1*q_2*intmaxWni_Win_h(i) ...
            +lif_1*intmaxWfn_Win_h(i) ...
            +lnf_2*intmaxWif_Wnf_Win_h(i) ...
            +lni_2*intmaxWii_Win_h(i)) + ...
            nu_hl*(1/(r+di_1*(1-q_2)))*(Win_l(i)-Win_h(i))+(1-nu_hl)*(1/(r+di_1*(1-q_2)))*0-Win_h(i);
  
        
        eWin_l(i)=(1/(r+di_1*(1-q_2)))*((1-exp(-theta*(wi_1(i)+b2_l)))/theta+gamma+di_1*(1-q_2)*Wnn_l ...
            +di_1*q_2*intmaxWni_Win_l(i) ...
            +lif_1*intmaxWfn_Win_l(i) ...
            +lnf_2*intmaxWif_Wnf_Win_l(i) ...
            +lni_2*intmaxWii_Win_l(i))+...
            nu_lh*(1/(r+di_1*(1-q_2)))*(Win_h(i)-Win_l(i))+(1-nu_lh)*(1/(r+di_1*(1-q_2)))*0-Win_l(i);
        
        
    end
    % Wnf_h(w2) and Wni_h(w2)
    for i=1:n
        
        
        eWnf_h(i)=(1/(r+df_2*(1-p_1)))*((1-exp(-theta*(wf_2(i)+b1_h)))/theta+a_h+df_2*(1-p_1)*Wnn_h ...
            +df_2*p_1*intmaxWin_Wnf_h(i) ...
            +lfi_2*intmaxWni_Wnf_h(i) ...
            +lnf_1*intmaxWff_Wfn_Wnf_h(i) ...
            +lni_1*intmaxWif_Wnf_h(i)) +...
            nu_hl*(1/(r+df_2*(1-p_1)))*(Wnf_l(i)-Wnf_h(i))+(1-nu_hl)*(1/(r+df_2*(1-p_1)))*0-Wnf_h(i);
      
        
        eWnf_l(i)=(1/(r+df_2*(1-p_1)))*((1-exp(-theta*(wf_2(i)+b1_l)))/theta+a_l+df_2*(1-p_1)*Wnn_l ...
            +df_2*p_1*intmaxWin_Wnf_l(i) ...
            +lfi_2*intmaxWni_Wnf_l(i) ...
            +lnf_1*intmaxWff_Wfn_Wnf_l(i) ...
            +lni_1*intmaxWif_Wnf_l(i))+...
            nu_lh*(1/(r+df_2*(1-p_1)))*(Wnf_h(i)-Wnf_l(i))+(1-nu_lh)*(1/(r+df_2*(1-p_1)))*0-Wnf_l(i);
        

        eWni_h(i)=(1/(r+di_2*(1-q_1)))*((1-exp(-theta*(wi_2(i)+b1_h)))/theta+gamma+di_2*(1-q_1)*Wnn_h ...
            +di_2*q_1*intmaxWin_Wni_h(i) ...
            +lif_2*intmaxWnf_Wni_h(i) ...
            +lnf_1*intmaxWfi_Wfn_Wni_h(i) ...
            +lni_1*intmaxWii_Wni_h(i))+...
            nu_hl*(1/(r+di_2*(1-q_1)))*(Wni_l(i)-Wni_h(i))+ (1-nu_hl)*(1/(r+di_2*(1-q_1)))*0-Wni_h(i);
        
        
        eWni_l(i)=(1/(r+di_2*(1-q_1)))*((1-exp(-theta*(wi_2(i)+b1_l)))/theta+gamma+di_2*(1-q_1)*Wnn_l ...
            +di_2*q_1*intmaxWin_Wni_l(i) ...
            +lif_2*intmaxWnf_Wni_l(i) ...
            +lnf_1*intmaxWfi_Wfn_Wni_l(i) ...
            +lni_1*intmaxWii_Wni_l(i))+...
            nu_lh*(1/(r+di_2*(1-q_1)))*(Wni_h(i)-Wni_l(i))+(1-nu_lh)*(1/(r+di_2*(1-q_1)))*0-Wni_l(i);
        
    end
    
    for i=1:n
        for j=1:n
            % Wff_h(w1,w2)
            
            eWff_h(i,j)=(1/(r+df_1+df_2))*((1-exp(-theta*(wf_1(i)+wf_2(j))))/theta+a_h+df_1*Wnf_h(j)+df_2*Wfn_h(i) ...
                +lfi_1*intmaxWif_Wff_h(i,j)...
                +lfi_2*intmaxWfi_Wff_h(i,j))+...
                nu_hl*(1/(r+df_1+df_2))*(Wff_l(i,j)-Wff_h(i,j))+(1-nu_hl)*(1/(r+df_1+df_2))*0-Wff_h(i,j);
            
            
            eWff_l(i,j)=(1/(r+df_1+df_2))*((1-exp(-theta*(wf_1(i)+wf_2(j))))/theta+a_l+df_1*Wnf_l(j)+df_2*Wfn_l(i) ...
                +lfi_1*intmaxWif_Wff_l(i,j)...
                +lfi_2*intmaxWfi_Wff_l(i,j))+...
                nu_lh*(1/(r+df_1+df_2))*(Wff_h(i,j)-Wff_l(i,j))+(1-nu_lh)*(1/(r+df_1+df_2))*0-Wff_l(i,j);
          
            
            % Wii_h(w1,w2)
            
            eWii_h(i,j)=(1/(r+di_1+di_2))*((1-exp(- theta*(wi_1(i)+wi_2(j))))/theta+gamma+di_1*Wni_h(j)+di_2*Win_h(i) ...
                +lif_1*intmaxWfi_Wfn_Wii_h(i,j) ...
                +lif_2*intmaxWif_Wnf_Wii_h(i,j))+...
                nu_hl*(1/(r+di_1+di_2))*(Wii_l(i,j)-Wii_h(i,j))+(1-nu_hl)*(1/(r+di_1+di_2))*0-Wii_h(i,j);

            
            eWii_l(i,j)=(1/(r+di_1+di_2))*((1-exp(- theta*(wi_1(i)+wi_2(j))))/theta+gamma+di_1*Wni_l(j)+di_2*Win_l(i) ...
                +lif_1*intmaxWfi_Wfn_Wii_l(i,j) ...
                +lif_2*intmaxWif_Wnf_Wii_l(i,j))+...
                nu_lh*(1/(r+di_1+di_2))*(Wii_h(i,j)-Wii_l(i,j))+(1-nu_lh)*(1/(r+di_1+di_2))*0-Wii_l(i,j);
           
            
            % Wfi_h(w1,w2)
            
            eWfi_h(i,j)=(1/(r+df_1+di_2))*((1-exp(-theta*(wf_1(i)+wi_2(j))))/theta+a_h+df_1*Wni_h(j)+di_2*Wfn_h(i) ...
                +lfi_1*intmaxWii_Wfi_h(i,j) ...
                +lif_2*intmaxWff_Wnf_Wfi_h(i,j))+...
                nu_hl*(1/(r+df_1+di_2))*(Wfi_l(i,j)-Wfi_h(i,j))+(1-nu_hl)*(1/(r+df_1+di_2))*0-Wfi_h(i,j);
            
            eWfi_l(i,j)=(1/(r+df_1+di_2))*((1-exp(-theta*(wf_1(i)+wi_2(j))))/theta+a_l+df_1*Wni_l(j)+di_2*Wfn_l(i) ...
                +lfi_1*intmaxWii_Wfi_l(i,j) ...
                +lif_2*intmaxWff_Wnf_Wfi_l(i,j))+...
                nu_lh*(1/(r+df_1+di_2))*(Wfi_h(i,j)-Wfi_l(i,j))+(1-nu_lh)*(1/(r+df_1+di_2))*0-Wfi_l(i,j);
           
            
            % Wif_h(w1,w2)
            
            eWif_h(i,j)=(1/(r+di_1+df_2))*((1-exp(-theta*(wi_1(i)+wf_2(j))))/theta+a_h+di_1*Wnf_h(j)+df_2*Win_h(i) ...
                +lif_1*intmaxWff_Wfn_Wif_h(i,j) ...
                +lfi_2*intmaxWii_Wif_h(i,j))+...
                nu_hl*(1/(r+di_1+df_2))*(Wif_l(i,j)-Wif_h(i,j))+(1-nu_hl)*(1/(r+di_1+df_2))*0-Wif_h(i,j);
           
            
            eWif_l(i,j)=(1/(r+di_1+df_2))*((1-exp(-theta*(wi_1(i)+wf_2(j))))/theta+a_l+di_1*Wnf_l(j)+df_2*Win_l(i) ...
                +lif_1*intmaxWff_Wfn_Wif_l(i,j) ...
                +lfi_2*intmaxWii_Wif_l(i,j))+...
                nu_lh*(1/(r+di_1+df_2))*(Wif_h(i,j)-Wif_l(i,j))+(1-nu_lh)*(1/(r+di_1+df_2))*0-Wif_l(i,j);
            
        end
    end
    % Wnn_h
    
    eWnn_h=(1/(r))*((1-exp(-theta*(b1_h+b2_h)))/theta+gamma ...
        +lnf_1*intmaxWfn_Wnn_h ...
        +lni_1*intmaxWin_Wnn_h ...
        +lnf_2*intmaxWnf_Wnn_h ...
        +lni_2*intmaxWni_Wnn_h)+...
        nu_hl*(1/(r))*(Wnn_l-Wnn_h)+(1-nu_hl)*(1/(r))*0-Wnn_h;
    
    
    eWnn_l=(1/(r))*((1-exp(-theta*(b1_l+b2_l)))/theta+gamma ...
        +lnf_1*intmaxWfn_Wnn_l ...
        +lni_1*intmaxWin_Wnn_l ...
        +lnf_2*intmaxWnf_Wnn_l ...
        +lni_2*intmaxWni_Wnn_l)+... 
        nu_lh*(1/(r))*(Wnn_h-Wnn_l)+(1-nu_lh)*(1/(r))*0-Wnn_l;
    
    % Update value functions
    
    Wfn_h=Wfn_h+omegaW*eWfn_h';
    Win_h=Win_h+omegaW*eWin_h';
    Wnf_h=Wnf_h+omegaW*eWnf_h';
    Wni_h=Wni_h+omegaW*eWni_h';
    Wff_h=Wff_h+omegaW*eWff_h;
    Wii_h=Wii_h+omegaW*eWii_h;
    Wfi_h=Wfi_h+omegaW*eWfi_h;
    Wif_h=Wif_h+omegaW*eWif_h;
    Wnn_h=Wni_h(1);eWnn_h=0;
    Wii_h(1,1)=Win_h(1);eWii_h(1,1)=0;
    
    Wfn_l=Wfn_l+omegaW*eWfn_l';
    Win_l=Win_l+omegaW*eWin_l';
    Wnf_l=Wnf_l+omegaW*eWnf_l';
    Wni_l=Wni_l+omegaW*eWni_l';
    Wff_l=Wff_l+omegaW*eWff_l;
    Wii_l=Wii_l+omegaW*eWii_l;
    Wfi_l=Wfi_l+omegaW*eWfi_l;
    Wif_l=Wif_l+omegaW*eWif_l;
    Wnn_l=Wni_l(1);eWnn_l=0;
    Wii_l(1,1)=Win_l(1);eWii_l(1,1)=0;
    
    critVF=[abs(eWfn_h'./Wfn_h);abs(eWin_h'./Win_h);abs(eWnf_h'./Wnf_h);abs(eWni_h'./Wni_h);max(abs(eWff_h./Wff_h))';max(abs(eWii_h./Wii_h))';max(abs(eWfi_h./Wfi_h))';...
        max(abs(eWif_h./Wif_h))';abs(eWnn_h/Wnn_h); abs(eWfn_l'./Wfn_l);abs(eWin_l'./Win_l);abs(eWnf_l'./Wnf_l)...
        ;abs(eWni_l'./Wni_l);max(abs(eWff_l./Wff_l))';max(abs(eWii_l./Wii_l))';max(abs(eWfi_l./Wfi_l))';max(abs(eWif_l./Wif_l))';abs(eWnn_l/Wnn_l)];
    
    iterVF=iterVF+1;
end
end
